Method of processing an input image by means of multi-resolution

ABSTRACT

The invention relates to a method of multi-resolution with gradient-adaptive filtering (MRGAF) of X-ray images in real time. For an image strip of 2K adjacent rows, a resolution into a Laplacian pyramid (L 0 , . . . L 3 ) and a Gaussian pyramid (G 0 , . . . G 3 ) is carried out up to the K-th stage. By limiting a processing operation to such a strip, it is possible to keep all relevant data ready in a local memory with rapid access (cache). A further acceleration compared to the conventional algorithm is achieved by calculating the gradient (D) from the Gaussian pyramid representations. By virtue of these and other optimization measures, it is possible to increase a multi-resolution with gradient-adaptive filtering to a processing rate of more than thirty (768 (E 564) images per second.

The invention relates to a method and a data processing unit forprocessing an input image, in particular for the multi-resolution andgradient-adaptive filtering of an X-ray image in real time.

The automated evaluation of images takes place in a large number ofdifferent fields of application. The processing of fluoroscopic X-rayimages considered in greater detail below should therefore be understoodto be merely one example. In order to minimize the amount of X-radiationto which a patient and the staff are subjected, X-ray images are takenwith as low a radiation dose as possible. However, there is a risk thatimportant image details will become lost in the image noise. In order toprevent this, attempts are made to suppress the noise by using spatialand temporal filters on the X-ray images or image sequences, withoutdestroying relevant image information in the process.

Within the context of such image processing, what is known asmulti-resolution of the input image is often carried out. The inputimage is in this case resolved into a sequence of detail images, wherethe detail images each contain image information from an associatedregion or strip at (spatial) frequencies. In addition, the detail imagesare adapted in terms of their resolution, i.e. the number of imagepoints for the representation of the content of the image, to theirrespective frequency range. By modifying detail images, it is possibleto have an influence on certain frequency ranges in a targeted manner.After the modification, the detail images can be put back together againto form an output image.

From WO 98/55916 A1 and EP 996 090 A2, in this respect an effectivemethod for the post-processing of X-ray images is known, in which amulti-resolution takes place and the detail images obtained thereby aremodified using filters, the coefficients of which have been adapted onthe basis of image gradients. The gradients in each case stem from thecoarser resolution stages of the multi-resolution. In this method, whichis referred to as MRGAF (Multi-Resolution Gradient Adaptive Filtering),low-pass filtering is carried out to a lesser extent perpendicular toimage structures such as lines or edges than in the direction of thestructures, so that a suppression of noise takes place that allowsinformation to be obtained. On account of the great calculation effortthat is required, however, the method can to date only be carried outoffline on stored images or image sequences.

In view of this background, it is an object of the present invention toprovide means for the more efficient processing of input images usingmulti-resolution, where preferably it should be possible to carry outimage analysis in real time.

This object is achieved by a method having the features as claimed inclaim 1, by a data processing unit having the features as claimed inclaim 8 and by an X-ray system having the features as claimed in claim10. Advantageous refinements are given in the dependent claims.

The method according to the invention is used for processing an inputimage comprising N rows of image points. Typically, the image points arearranged in a rectangular grid having columns perpendicular to the rows,although other arrangements having a row structure, such as e.g.hexagonal grids, are also possible. The input image may in particular bea digitized fluoroscopic X-ray image, although the method is notrestricted to this and can be advantageously used in all comparableapplication instances in which a multi-resolution of an image takesplace. The method comprises the following steps:

a) An image strip which comprises n<N adjacent rows of the input imageis resolved into a sequence of K detail images, where the detail imagesin each case contain just a partial range of the spatial frequencies ofthe image strip. A multi-resolution is thus carried out using thestrip-like section of the overall input image.

b) At least one of the detail images obtained in step a) is modified,for example using a predefined filter or a filter that is calculatedfrom the image strip. Preferably, all the information that is requiredfor the modification is available in the image strip.

c) An output image strip is reconstructed from the detail images or themodified detail images (if the latter exist).

d) The above steps a), b) and c) are repeated for other image strips ofthe input image, that is to say they are carried out in an analogousmanner with calculation of a corresponding output image strip. Othervalues for the strip width n and/or the resolution depth K may also beassumed where appropriate.

e) An output image is reconstructed from the calculated output imagestrips.

As a result, in the above method, an output image is therefore generatedfrom an input image (said output image being of the same size or of adifferent size), where modifications of the desired type have beencarried out in some or all spatial frequency ranges of the input image.Compared with conventional multi-resolution with a modification ofdetail images, the difference with this method is that themulti-resolution takes place in sections on image strips of in each casen rows. Each image strip is in this case resolved to the stage K andthen synthesized again to give an output image strip. The advantage ofthis procedure is that it is particularly suitable for efficientimplementation on a data processing system, since the memory requirementfor the processing of an image strip is accordingly smaller than for theprocessing of the full image, so that the method can be carried outusing a working memory with rapid access. As a result, a gain in speedcan be achieved that is so great that in many cases for the first timeit is practical for the multi-resolution to be carried out in real time.

According to one specific refinement of the method, in themulti-resolution of step a), each image strip is resolved into aLaplacian pyramid and a Gaussian pyramid with in each case K stages. Instage j of a Gaussian pyramid, the stage input image is the output imageof the preceding stage (j−1), and the output image (hereinafter referredto as “Gaussian pyramid representation of stage j”) is the stage inputimage that has been modified by low-pass filtering and subsequentresolution reduction. The output image of the Laplacian pyramid at stagej (hereinafter referred to as “Laplacian pyramid representation of stagej”) is obtained by subtracting the Gaussian pyramid representation ofthe same stage j, the resolution of which has been increased again andwhich has been low-pass-filtered, from the Gaussian pyramidrepresentation of the preceding stage (j−1). The resolution of an inputimage into a Laplacian pyramid or Gaussian pyramid is frequently used inmedical image processing and is particularly suitable for use on imagestrips.

Preferably, in steps a) to d), the image strip subjected tomulti-resolution is in each case 2^(K) rows wide, where K is the numberof resolution stages of the multi-resolution. Image strips having awidth of 2^(K) have the minimum width necessary for resolution into aLaplacian pyramid or Gaussian pyramid to the stage K, if at each stageof the resolution a reduction of the rows and columns by in each case afactor of 2 takes place. The detail image of the coarsest stage has theminimum width of one row for such an image strip. Furthermore, the imagestrips are optionally offset by in each case (2^(K)−1) rows with respectto one another, or in other words overlap one another by in each caseone row. Such an overlapping, which preferably is also present on allresolution stages of the image strips, provides the necessaryinformation for filter operations taking place at the edge of the new,non-overlapping region. Depending on the width of the filter used, theremay also be instances of overlapping of more than one row wide betweenthe image strips.

The type of modifications carried out with the detail images may differdepending on the application instance. Preferably, the modification of adetail image of the resolution stage j<K is in the use of a filter,where the coefficients of this filter depend on at least one gradientcalculated from the image strip. Since gradients of the image reflectthe position of local structures in the image, they can be used todefine anisotropic filters, the use of which leaves the structuresunchanged or even amplifies them, and suppresses any noise along thestructures.

Preferably, the above method is combined with a resolution into aGaussian pyramid and a Laplacian pyramid, and the gradient is calculatedfrom the Gaussian pyramid representation of the resolution stage j andis used for the filtering of the Laplacian pyramid representation of thesame stage j. This has the advantage that all information required forthe modification can be obtained from the data of the resolution stagej, so that the modification can be carried out directly during thecalculation of this stage.

According to a special design of the above gradient-adaptive filteringof detail images, the filter coefficients α(Δ{right arrow over(x)},{right arrow over (x)}) are calculated from the coefficientsβ(Δ{right arrow over (x)}) of a predefined filter, such as for example abinomial filter, where {right arrow over (x)} is the image pointprocessed by the filter and Δ{right arrow over (x)} is the position ofthe respective coefficient in relation to the center of the filter, andwhere the following formula applies:α(Δ{right arrow over (x)},{right arrow over (x)})=β(Δ{right arrow over(x)})[r({right arrow over (g)}({right arrow over (x)}), Δ{right arrowover (x)})]²  (1)

Herein, {right arrow over (g)}({right arrow over (x)}) is the gradientat the image position {right arrow over (x)} and 0≦r≦1. Where theweighting function r({right arrow over (g)},{right arrow over (x)})<1,the corresponding filter coefficients β are decreased and thecontribution thereof to the result of filtering is reduced. In this way,a noise contribution is suppressed at the corresponding positions of animage.

The weighting function r is preferably defined as follows:$\begin{matrix}{{r\left( {\overset{->}{g},{\Delta\quad\overset{->}{x}}} \right)} = \left( \frac{1}{1 + {{c\left\lbrack \overset{->}{g} \right\rbrack}\left( {{\overset{->}{g} \cdot \Delta}\quad\overset{->}{x}} \right)^{2}}} \right)} & (2)\end{matrix}$

-   -   where c[{right arrow over (g)}] is a factor that preferably        depends on the gradient field {right arrow over (g)} and the        variance thereof. The above definition of the factor r has the        desired property that r=1 in directions Δ{right arrow over (x)}        perpendicular to the gradient {right arrow over (g)}({right        arrow over (x)}), and that r is minimal in directions Δ{right        arrow over (x)} parallel to {right arrow over (g)}({right arrow        over (x)}). The definitions for the calculation of α and r are        considerably easier in terms of their calculation effort than        the definitions given in WO 98/55916 A1, with the results being        approximately identical.

The invention furthermore relates to a data processing unit forprocessing a digital input image comprising N rows of image points,which data processing unit contains a system memory and a cache memory.The data processing unit is intended to carry out the followingprocessing steps:

a) resolution of an image strip comprising n<N adjacent rows of theinput image into a sequence of K detail images, which in each casecontain just some of the spatial frequencies of the input image;

b) modification of at least one of the detail images;

c) reconstruction of an output image strip from the—possiblymodified—detail images;

d) repetition of steps a), b) and c) for other image strips of the inputimage;

e) reconstruction of an output image from the calculated output imagestrips;

wherein during steps a)-c) all processed data (data of the image strips,data of the associated detail images of the multi-resolution of theimage strips) is in each case located in the cache memory.

Using such a data processing unit, the above-described method can becarried out very efficiently and quickly, since all the necessary datacan be accommodated in the cache memory and thus can be accessedrapidly. By contrast, in conventional multi-resolution, in each case thefull image is analyzed, as a result of which use must be made of thesystem memory (working memory, hard disk, etc.) for the storage of theintermediate results. A large part of the calculation time is thus takenup with the reading and writing of the data to and from the systemmemory. Because these time-consuming operations are omitted, it ispossible using the above data processing unit to carry out the imageprocessing even in real time.

Preferably, the data processing unit is equipped with parallelprocessors and/or vector processors. In this case, the necessaryoperations can be speeded up even further by means of parallelization.

Furthermore, the data processing unit is preferably designed such thatit can also carry out the variants of the method explained above.

The invention further relates to an X-ray system comprising:

an X-ray source;

an X-ray detector;

a data processing unit, coupled to the X-ray detector, for processingthe X-ray input images generated by the X-ray detector, where the dataprocessing unit is designed in the above-described manner.

The advantage of such an X-ray system is that effective image processingcan be carried out in real time, i.e. during a medical intervention,said image processing suppressing noise without impairing structures ofinterest. In particular, an MRGAF method can be carried out in realtime. On account of the suppression of noise, which allows informationto be obtained, it is possible to take X-ray photographs with acorrespondingly low dose of radiation, and thus to minimize the amountof radiation to which the patient and the staff are subjected.

The invention will be further described with reference to examples ofembodiments shown in the drawings to which, however, the invention isnot restricted.

FIG. 1 shows the sequence of an MRGAF algorithm according to the priorart;

FIG. 2 shows the use of variables in the low-pass filtering andresolution reduction during the generation of a Gaussian pyramidrepresentation of the next-higher resolution stage;

FIG. 3 shows the calculation of the Laplacian pyramid representation andof the gradient fields in the x- and y-direction from two successiveGaussian pyramid representations;

FIG. 4 shows the position of the image points in various resolutionstages;

FIG. 5 shows the position of the image points in various compositionstages;

FIG. 6 shows the sequence of an MRGAF algorithm according to theinvention.

The MRGAF algorithm shown schematically in FIG. 1 is described in detailin EP 996 090 A2 and WO 98/55916 A1 and shall therefore only bedescribed by way of an overview below. The aim of the MRGAF algorithm isto significantly reduce the noise in an input image I while at the sametime maintaining the image details and the image sharpness. The basicidea of the algorithm consists in a multi-resolution and an anisotropiclow-pass filtering of the resulting detail images as a function of thelocal image gradient.

In the example shown in FIG. 1, resolution of the input image I, whichcomprises 512×512 image points (pixels), takes place in K=4 resolutionstages. At each resolution stage j=0, 1, 2, 3, what is known as aLaplacian pyramid representation Λ_(j) and a Gaussian pyramidrepresentation Γ_(j) is generated as detail image. The stage inputrepresentation is in each case the Gaussian pyramid representationΓ_(j−1) of the preceding stage (j−1) or the original inputrepresentation L The Gaussian pyramid representation Γ_(j) is generatedby using a reduction operation R on the respective stage inputrepresentation, where a “reduction” means a low-pass filtering(smoothing) and subsequent resolution reduction (subsampling) by thefactor 2, which leads to an image of half the size. The Laplacianpyramid representations Λ_(j) are defined as the difference between thestage input representation and the copy thereof after passing throughthe reduction R and expansion E blocks. The “expansion” E here includesa resolution increase by the factor 2 (by inserting zeros) and asubsequent low-pass filtering (interpolation). In this case, 3×3binomial filters are used for the low-pass filtering operations in thereduction R and the expansion E. The Laplacian pyramid representationsΛ_(j) accordingly contain the high-pass fraction and the Gaussianpyramid representations Γ_(j) contain the associated low-pass fractionof the resolution stage j (cf. B. Jähne, Digitale Bild-verarbeitung[Digital Image Processing], 5th edition, Springer Verlag BerlinHeidelberg, 2002, Section 11.4, 5.3).

In preparation for the gradient filtering, by means of simple differenceformation between adjacent pixels the gradients Δ are furthermorecalculated from the Laplacian pyramid representations Λ_(j). Therespective difference in this case belongs to a location in the centerbetween the pixels used for difference formation. Furthermore, althoughthe gradient is calculated at the resolution stage j, it is used forfiltering at the preceding, finer resolution stage (j−1). For thesereasons, the gradients have to be suitably interpolated. Finally, theresult is again divided by the factor 2, in order to compensate for thefiner sampling. Since the modulus of the band-pass image does not onlyassume maximum values in the event of discontinuities in the originalimage but also in the vicinity thereof, the gradients of the coarserresolution stages j′>j are expanded in the block E and added to thegradient of the resolution stage j.

With the exception of the coarsest resolution stage, all Laplacianpyramid representations Λ₀ to Λ₂ are filtered using a filter GAF, whichreacts adaptively to the gradients calculated as described. The startingpoint of the filter synthesis is in this case a 3×3 binomial filter, thefilter coefficients β(Δ{right arrow over (x)}) of which are maintainedalong the main directions of the representation structures, while thecoefficients in the direction of the gradients to these structures arereduced according to the following formula:α(Δ{right arrow over (x)},{right arrow over (x)})=β(Δ{right arrow over(x)})r({right arrow over (g)}({right arrow over (x)}),Δ{right arrow over(x)})r({right arrow over (g)}({right arrow over (x)}+Δ{right arrow over(x)}),Δ{right arrow over (x)})  (3)

-   -   where α(Δ{right arrow over (x)},{right arrow over (x)}) is the        new coefficient of the filter, {right arrow over (x)} is the        image position to be filtered, Δ{right arrow over (x)} is the        vector pointing from the center of the filter core to the        coefficient in question, β(Δ{right arrow over (x)}) is the        original filter coefficient, r is a weighting function and        {right arrow over (g)}({right arrow over (x)}) is the gradient        at the image point {right arrow over (x)}. The weighting        function r drops exponentially with the scalar product of the        gradient and of the coefficient direction Δ{right arrow over        (x)} as shown in $\begin{matrix}        {{r\left( {\overset{->}{g},{\Delta\quad\overset{->}{x}}} \right)} = {\exp\left( {- \frac{\left( {{\overset{->}{g} \cdot \Delta}\quad\overset{->}{x}} \right)^{2}}{c + {t \cdot {{Var}\left( \overset{->}{g} \right)}} + {L{\overset{->}{g}}^{2}}}} \right)}} & (4)        \end{matrix}$    -   where c, t and L are parameters that can be defined by the user        and Var({right arrow over (g)}({right arrow over (x)})) is the        estimated variance of the noise of the gradient field.

As can be seen in FIG. 1, the calculation of the gradient-adaptivefilter GAF presupposes the already processed and reconstructed Gaussianpyramid with the representations Γ_(j).

The right-hand part of the diagram in FIG. 1 reflects the synthesis ofan output image A from the detail images Λ_(j) (which are unmodified orhave been modified by a filtering GAF) by means of successive additionand expansion E. If no filtering of the detail images were to have takenplace, the output image A would be identical to the input image I.

A disadvantage of the described MRGAF algorithm is that to date it canonly be carried out offline on stored images or image sequences onaccount of the high calculation effort. Because of the significant imageimprovement that can be achieved with this algorithm, however, it wouldbe desirable to be able to carry it out also in real time, for exampleduring an ongoing medical intervention. This aim is achieved in themanner described below using various optimizations, but particularly byan approach for processing what are known as “super-rows”. Thisprocessing principle cannot only be used in the case of the MRGAFalgorithm considered here by way of example; rather it can be used inprinciple in all types of multi-resolution and also with othercomparable algorithms, such as for example a “sub-band coding”.

The above-described original MRGAF algorithm processes the data in alevel by level manner. First, the input image I is low-pass-filtered.Since the images are typically too large to pass into a (buffer) memorywith rapid access by the processor (cache), some of the input data andsome of the processed data must be read from or written to the main orsystem memory (working memory RAM and/or bulk memory, such as e.g. harddisk). However, this is disadvantageous for two reasons: firstly, theaccesses to the system memory are relatively slow and, secondly, memoryhardware does not progress as rapidly in technological terms as dataprocessing hardware with regard to the increase in speed. Therefore, ithas been sought to change the MRGAF algorithm in such a way that thenumber of memory accesses is reduced.

In order to reduce the number of read/write operations on the systemmemory, the calculations of the Gaussian pyramid representations and ofthe Laplacian pyramid representations and also of the gradient fields ateach resolution stage j are combined with one another, so that thesevalues are calculated locally in a single pass over the stage inputimage. For this purpose, the low-pass-filtered and resolution-reducedvalue of the next-coarser Gaussian pyramid representation Γ_(j+1) mustfirst be calculated. The fastest way to do this is shown schematicallyin FIG. 2. Instead of using a 3×3 binomial low-pass filter (1,2,1;2,4,2; 1,2,1)· 1/16, multiplication and addition operations can be savedby using the one-dimensional low-pass filter (1,2,1)·¼ successively inthe x- and y-direction (that is to say in the row and column direction)and by buffering the intermediate values b_(i). Specifically, thealgorithm shown in FIG. 2 for calculating a value from Γ_(j+1) (point inFIG. 2) proceeds as follows: Let x_(0,0) x_(0,1) x_(0,2) x_(1,0) x_(1,1)x_(1,2) x_(2,0) x_(2,1) x_(2,2) be a 3 × 3 proximity around the currentposition x_(1,1) and b₁ = x_(0,0) + x_(0,1) + x_(0,1) + x_(0,2) be avalue buffer-stored from the preceding row.

The following steps are then carried out, where v₁, v₂ are temporaryvariables and b₀, b₁, b₂, . . . are buffer variables:

-   -   1. v_(i)=x_(1,0)+x_(1,1)+x_(1,1)+x_(1,2)    -   2. v₂=x_(2,0)+x_(2,1)+x_(2,1)+x_(2,2)    -   3. result= 1/16*(b₁+v₁+v₁+v₂)→enter in Γ_(j+1)    -   4. b₁=v₂    -   (etc.)

Furthermore, as shown in FIG. 3, the resulting value of the Gaussianpyramid representation Γ_(j+1) is used directly to calculate in eachcase four values in the Laplacian pyramid representation Λ_(j) and inthe gradient fields Δ_(x) in the x-direction and Δ_(y) in they-direction (cf. marked points in FIG. 3). The Laplacian values areobtained here from the subtraction

of the current value of the Gaussian pyramid representation Γ_(j+1) andof the values interpolated herewith in relation to the alreadycalculated neighbors in Γ_(j+1) to the left of and above the currentposition

of the corresponding values of the Gaussian pyramid representationΓ_(j).

The gradient values are the difference from the already calculatedvalues of the Gaussian pyramid representation Γ_(j+1) to the left of andabove the current position (short dash in FIG. 3) and the interpolatedvalues using the already calculated values in the gradient fields. Byvirtue of the resolution increase that is carried out at the same time,it is possible to set the calculated differences at the correct“intermediate positions”.

Using the described memory-optimized calculation of the data requiredfor the GAF routine, the MRGAF algorithm can already be carried outaround 15% quicker. A further significant increase in efficiency isachieved by the innovation that the overall resolution is carried outwith the smallest possible amount of data, so that the data required inthis case can be buffered in a memory with rapid access (cache). In thiscase, the read/write operations for the input image and the output imageare the only accesses to the slower system memory that are stillrequired.

Since a read/write command at a memory address always leads to thereading/writing of the entire subsequent data block from or to thecache, the processing of complete rows is retained. However, it is notthe entire input image I which is processed in one go, but rather onlyas few rows as possible. This means that, as shown in FIG. 4, theGaussian pyramid representation Γ₃ of the coarsest resolution stage onlycomprises a single row together with the preceding row, which isrequired for interpolation and difference calculation. Therefore, onaccount of the pyramid structure, a “super-row” of 2^(K) rows must beprocessed at the same time, where K is the maximum resolution stage(e.g. K=3 in FIGS. 4, 5).

FIG. 3 can be seen as a representation of the calculation of theLaplacian pyramid block and of the gradient blocks at the coarsestresolution stage. The blocks comprise two rows plus an additionalpreceding row for the interpolation. As can be seen in FIG. 3, adisplacement takes place on account of the reduction, the re-expansionand the interpolation. If the relative rows 0 and 1 are given as input,the gradient-adaptive filtering GAF can be carried out only on the rows−1 and −2 of the Laplacian pyramid block. The reason for this is theposition of the resulting data of the y-gradient and the fact that thefiltering with a 3×3 GAF filter core requires a pixel of additional dataat each side of the filter position. This additional data is the rows 0and −3 (not shown in FIG. 3) and the first and last columns of theLaplacian pyramid block.

FIG. 4 shows how the displacement effect responds at the otherresolution stages. It can be seen that the filtered area (dark gray)always lies two rows above the current position and the Laplacianpyramid block Λ_(j) (light gray) lies one row above the currentposition, where the latter is expanded at the top by two preceding rowsin order to permit a 3×3 filtering operation.

In the last step of the method, the image block is reconstructed fromthe filtered Laplacian pyramid representations. As shown in FIG. 5, thedisplacement of the filtered data by two rows is summed during thesynthesis step, and this results in a relatively large displacement ofthe reconstructed image block. However, this does not mean that the datahave been rewritten at the wrong points; all values pass back to wherethey came from. In reality, it is not a displacement in terms oflocation, but rather in terms of time on account of the causalitycondition that means interpolation can only take place with alreadycalculated values. The light gray areas in FIG. 5 therefore distinguishpreviously filtered data that has to last until synthesis. In thisrespect, however, a considerable amount of buffer data can be saved bysimply swapping the steps of filtering and synthesis: firstly, the onlythree rows that are still required for synthesis are filtered (two rowsfrom resolution stage 2 and one row from stage 1—dark gray in FIG. 5).Then the synthesis is carried out, so that the data in thereconstruction buffers can be overwritten with the remaining filtervalues (dark gray). As a result of this swapping, the reconstructionbuffer of stage 0 for example does not need to keep nine additionalpreceding rows ready, but rather just one. This number would be 3, 5, 7,. . . if more resolution stages were used.

The following calculation estimates the overall memory requirement forbuffering data in the above method. It is assumed here that the imagewidth is 512 pixels, a three-stage pyramid is used, and 4 byte floatingdecimal values are used for all calculations. $\begin{matrix}{{Gaussian}\quad{{pyramid}:}} & {4*\left\lbrack {{512*8} + {256*\left( {4 + 1} \right)} + {128*\left( {2 + 1} \right)} + {64*\left( {1 + 1} \right)}} \right\rbrack} & {= {23\quad 552\quad{bytes}}} \\{{Laplacian}\quad{{pyramid}:}} & {4*\left\lbrack {{512*\left( {8 + 2} \right)} + {256*\left( {4 + 2} \right)} + {128*\left( {2 + 2} \right)}} \right\rbrack} & {= {28\quad 672\quad{bytes}}} \\{{Gradient}\quad{{fields}:}} & {2*4*\left\lbrack {{512*\left( {8 + 1} \right)} + {256*\left( {4 + 1} \right)} + {128*\left( {2 + 1} \right)}} \right\rbrack} & {= {50\quad 176\quad{bytes}}} \\{{Synthesis}\quad{{buffer}:}} & {4*\left\lbrack {{512*\left( {8 + 1 + 1} \right)} + {256*\left( {4 + 1} \right)} + {128*\left( {2 + 1} \right)}} \right\rbrack} & {= {27\quad 136\quad{bytes}}} \\\quad & \quad & {{\Sigma 129}\quad 536\quad{bytes}}\end{matrix}\quad$

The calculation shows that all calculations can be carried out usingdata in the second-level cache if the latter has a typical size of 256kB=262 144 bytes. There is even still space for a further 19 rows of theoriginal image if the result is to be rewritten hereto:

-   -   4*19*512=38912 bytes

For a comparison of the above methods, the original version of the MRGAFalgorithm can be sketched as follows:

For all levels of the pyramid:

-   1. Reduction of the stage input image in the x-direction→buffer-   2. Reduction of the buffer in the y-direction→Gaussian pyramid    representation-   3. Expansion of the Gaussian pyramid representation in the    y-direction→buffer-   4. Expansion of the buffer in the x-direction→second buffer-   5. Subtraction of the second buffer from the stage input    image→Laplacian pyramid representation

By contrast hereto, in the case of the new algorithm, everything takesplace using data from the second-level cache during just a single passof the image. The pyramid resolution, the gradient calculation, theadaptive filtering and the image synthesis are carried out on imagestrips that are as narrow as possible.

The size of the temporary buffer memory is derived from the requirementthat the coarsest stage of the Gaussian pyramid comprises only one rowtogether with the preceding row for the interpolation. The situation iscomplicated by the fact that, on account of all re-expansions withinterpolations, which are only possible using already processed rows,the filtering can be carried out no further than up to two rows abovethe lower limit of the read image block (cf. FIG. 4: the y-gradientscannot be calculated for the last two rows). For the coarsest stage ofthe Laplacian pyramid with a height of two pixels, this means that it isonly possible to filter the preceding block. During the reconstruction,this displacement grows from stage to stage. This means that a largeamount of data (at least the size of one block) needs to be kept readyand displaced in the buffers in each step. However, by swapping the“filtering” and the “reconstruction”, the number of copying operationscan fortunately be reduced. The algorithm is shown schematically in FIG.6, said algorithm comprising the following steps:

For all image strips of the

-   -   size 2ˆ(number of resolution stages)×(image width):    -   1. For all resolution stages:

Passes of the stage image strips in the x- and y-direction in the 2ndsteps, where the following is carried out at each position:

low-pass filtering in the x- and y-direction→one pixel of the Gaussianpyramid representation

expansion of the written pixel back to four pixels (using its neighborfor the interpolation) and direct subtraction thereof from the pixels ofthe stage input representation→4 pixels of the Laplacian pyramidrepresentation

calculation of the difference between the calculated pixel of theGaussian pyramid representation and its neighbor, interpolation of theresults with previously calculated gradients→4 pixels in the x- andy-gradient fields

2. Filtering of the last two rows of the coarsest stage of the Laplacianpyramid and of the first row from the second-coarsest stage (these rowsare still required for the reconstruction, the others are alreadyavailable)→reconstruction buffer

3. Reconstruction of an image strip from the reconstruction pyramidbuffer

4. For all stages other than the coarsest:

copying of the data required in the next step from bottom to top in thereconstruction buffer

filtering of the current Laplacian pyramid strip→reconstruction buffer.

In the text which follows, a series of refinements of the algorithm aredescribed, which contribute to further acceleration.

From a comparison of FIG. 1 with FIG. 6, which shows the MRGAF algorithmaccording to the invention, an important difference can be seen in thefact that in the new algorithm the gradient fields are calculated ateach stage directly from the low-pass-filtered stage inputrepresentations. Since the algorithm no longer needs to wait until thenext-coarser pyramid stage is filtered, the filtering is now carried outin parallel at all resolution levels.

As shown in equation (4), in the original MRGAF algorithm an exponentialfunction is calculated for each filter coefficient at each filterposition and at each resolution stage. In this respect, it is proposedto replace the function exp(−x) with the approximation 1/(1+x), whichprovides a similar profile for a considerably lower calculation effort.In this way, one arrives at equation (2) for the new filtercoefficients.

In the original MRGAF algorithm, as shown in equation (3) each filtercoefficient is weighted with two factors r, of which one is calculatedwith the gradients at the current image position {right arrow over (x)}and the other is calculated with the gradients at the coefficientposition {right arrow over (x)}+Δ{right arrow over (x)}. In the case ofa 3×3 filter core, therefore, nine different gradient factors have to becalculated for each filtered pixel. By virtue of the gradientcalculation at all filter positions, the treatment of curved lines oughtto be improved. However, when using 3×3 filters, this procedure provesto be unnecessary, since adjacent gradients interpolated from thecoarser pyramid stage are very similar to one another. Instead ofequation (3), therefore, the simplified formula as shown in equation (1)is proposed. The calculation of the filter routine is considerablysimplified on account of this, since opposite filter coefficients nowhave the same value and equation (2) needs to be calculated only once,instead of nine times.

The variance of the noise of the gradient field in the denominator ofequation (2), Var({right arrow over (g)}({right arrow over (x)})), isapproximately replaced by the variance of the noise of the correspondingpixel of the coarser pyramid stage. The quality of the filter result isnot impaired by this.

From FIG. 6 it can be seen that, when using the Gaussian pyramidrepresentations for gradient calculation, the coarsest pyramid stage(with Γ₃, Λ₃) is superfluous. The calculation of this stage cantherefore be dispensed with. A three-stage filtering operation thusleads to the same result as was previously achieved with a four-stagefiltering operation.

On a dual Xeon Pentium 4 with 1.7 GHz and with compilation using anIntel C++ compiler, in the most favorable case the implementation of aprogram taking account of the simplifications discussed above led to arun time of 0.0229 sec for one image, corresponding to 43.6 images(512×512) per second (i.e. more than thirty (768×564) images/s). Themethod has thus reached a point such that it is suitable for real-timeapplications.

1. A method of processing an input image (I) comprising N rows of imagepoints, wherein a) an image strip comprising n<N adjacent rows of theinput image is resolved into a sequence of K detail images (Λ₀, . . .Λ₃; Γ₀, . . . Γ₃), which in each case contain just some of the spatialfrequencies of the input image; b) at least one of the detail images(Λ₀, . . . Λ₂) is modified; c) an output image strip is reconstructedfrom the—possibly modified—detail images; d) steps a), b) and c) arerepeated for other image strips of the input image; e) an output image(A) is reconstructed from the calculated output image strips.
 2. Amethod as claimed in claim 1, characterized in that each image strip isresolved into a Laplacian pyramid and a Gaussian pyramid with K stages.3. A method as claimed in claim 1, characterized in that the imagestrips each have a width of 2^(K) rows.
 4. A method as claimed in claim1, characterized in that the modification of a detail image (Λ_(j)) ofthe resolution stage j<K is effected using a filter (GAF), thecoefficients of which depend on at least one gradient calculated fromthe image strip.
 5. A method as claimed in claims 2 and 4, characterizedin that the gradient is calculated from the Gaussian pyramidrepresentation (Γ_(j)) of the j-th resolution stage.
 6. A method asclaimed in claim 4, characterized in that the filter coefficientsα(Δ{right arrow over (x)},{right arrow over (x)}) are calculated fromthe coefficients β(Δ{right arrow over (x)}) of a predefined filter,according to the formulaα(Δ{right arrow over (x)},{right arrow over (x)})=β(Δ{right arrow over(x)})[r({right arrow over (g)}({right arrow over (x)}),Δ{right arrowover (x)})]² where {right arrow over (x)} is the image point processedby the filter operation, Δ{right arrow over (x)} is the position of thecoefficient relative to the center of the filter, {right arrow over(g)}({right arrow over (x)}) is the gradient at the image position{right arrow over (x)} and 0≦r≦1.
 7. A method as claimed in claim 6,characterized in that${{r\left( {\overset{->}{g},{\Delta\quad\overset{->}{x}}} \right)} = \left( \frac{1}{1 + {{c\left\lbrack \overset{->}{g} \right\rbrack}\left( {{\overset{->}{g} \cdot \Delta}\quad\overset{->}{x}} \right)^{2}}} \right)},$where c[{right arrow over (g)}] is a positive factor that is preferablydependent on the gradient field and its variance.
 8. A data processingunit for processing a digital input image (I) comprising N rows of imagepoints, which data processing unit contains a system memory and a cachememory and is intended to carry out the following processing steps: a)resolution of an image strip comprising n<N adjacent rows of the inputimage into a sequence of K detail images (Λ₀, . . . Λ₃; Γ₀, . . . Γ₃),which in each case contain just some of the spatial frequencies of theinput image; b) modification of at least one of the detail images (Λ₀, .. . Λ₂); c) reconstruction of an output image strip from the—possiblymodified—detail images; d) repetition of steps a), b) and c) for otherimage strips of the input image; e) reconstruction of an output image(A) from the calculated output image strips; wherein during steps a)-c)all processed data is in each case located in the cache memory.
 9. Adata processing unit as claimed in claim 8, characterized in that itcontains parallel processors and/or vector processors.
 10. An X-raysystem comprising an X-ray source; an X-ray detector; a data processingunit as claimed in claim 8, coupled to the X-ray detector, forprocessing the X-ray input images transmitted by the X-ray detector.